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In this paper we propose a novel method to study critical systems numerically by a combined 
collective-mode algorithm and Renormalization Group on the lattice. This method is an improved 
version of MCRG in the sense that it has all the advantages of cluster algorithms. As an application 
we considered the 2D Ising model and studied wether scale invariance or universality are possible 
underlying mechanisms responsible for the approximate "universal fluctuations" close to a so-called 
bulk temperature T*(L). "Universal fluctuations" was first proposed in [lj and stated that the 
probability density function of a global quantity for very dissimilar systems, like a confined turbulent 
flow and a 2D magnetic system, properly normalized to the first two moments, becomes similar to 
the "universal distribution", originally obtained for the magnetization in the 2D XY model in the 
low temperature region. The results for the critical exponents and the renormalization group flow 
of the probability density function are very accurate and show no evidence to support that the 
approximate common shape of the PDF should be related to both scale invariance or universal 
behavior. 
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I. INTRODUCTION 

Critical phenomena are present in a large number 
of quite different physical systems: super-fluid Helium 
three, low temperature super-conductors, ferromagnetic- 
paramagnetic systems, turbulent fluids, plasmas, poly- 
mers, among many others. Nevertheless, an important 
common feature of these systems is the scale indepen- 
dent fluctuations at the critical temperature: although 
the underlying inter-molecular forces, responsible for the 
existence of phase transitions, have a well-defined length 
scale, the structures they give rise do not. This leads, 
very close to the critical temperature (scaling region) , to 
a power-law behavior of the physical quantities, which is 
a fundamental feature of universality [4] . 

The main challenge of the theory of critical phenom- 
ena is to explain how dissimilar systems exhibit the same 
critical behavior. Renormalization Group (RG) Theory 
developed by Wilson and Kogut [l| provides a consistent 
framework to understand the existence of equivalence 
classes of critical systems. On the lattice, a very cun- 
ning method of applying a RG analysis to Monte Carlo 
simulations of general systems was first proposed in (4)]. 

Nevertheless, the numerical simulation of critical sys- 
tems has a serious limitation due to the critical slowing 
down effect. Indeed, as a critical system approaches the 
critical temperature, the decorrelation time diverges with 
the power of the correlation length of the system £ to the 
dynamical critical exponent z: r ~ £ z , where £ is ap- 
proximately 2 for local-flip algorithms like the Metropolis 
algorithm. In order to beat or at least reduce this effect, 
a cluster algorithm was first developed in Q . 
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In this paper, we develop a self-consistent method 
along the line proposed by Swendsen in which will 
be explained in section II, to perform a Lattice Renor- 
malization Group (RG) analysis of the probability den- 
sity function (PDF) as a function of the magnetization, 
in the 2D Ising model. His proposal amounts to compute 
by a direct Monte Carlo simulation of the fundamental 
Hamiltonian, a sequence of approximations to the lin- 
earized RG matrix T a p* . From its eigenvalues the crit- 
ical exponents can be obtained in a direct way. In or- 
der to reduce critical slowing down, we used instead of 
a local (Metropolis-type) algorithm the collective-mode 
algorithm developed in [6] to simulate the fundamental 
Hamiltonian. We call this method Cluster Algorithm 
Renormalization Group (CARG) [7[. 

Our physical motivation is to study, in the context of 
the 2D Ising model, whether scale invariance and univer- 
sality are the underlying mechanism which could give rise 
to the approximate "universal curve" . This phenomenon 
link a large class of dissimilar systems defined on differ- 
ent dimensions and including non-equilibrium systems, 
in which the PDF of a global quantity -like the power 
consumption in a confined turbulent system or the mag- 
netization in a finite ferromagnetic system- properly nor- 
malized to the first two moments, is described by a sin- 
gle curve Q. This claim is not free of controversy, and 
in the last years it has been the central task of several 
publications [9l]-[l6j|. This "universal curve" was shown 
to correspond to the PDF of the 2D-XY model in the 
zero temperature limit (T3 |. 

The paper is organized as follows: in section II the 
CARG method is explained for a general critical system. 
In the third section we explain the concept of approxi- 
mate "universal fluctuations" of the PDF of the magne- 
tization for the 2D Ising model. In section IV the criti- 
cal exponents v and r\ are computed using fundamental 
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lattices of lattice sizes L = 45, L = 64 and L = 108, per- 
forming two and three RG steps, to check the accuracy of 
the method. It is shown that relative errors of the criti- 
cal exponents increase monotonically with the departure 
from the critical temperature. Finally, the RG flow of 
the PDF itself is computed and it is shown that when 
conveniently normalized, it is invariant. The underlying 
reason for this invariance is discussed in connection with 
the so-called "generalized universality" conjecture. The 
results are discussed and conclusions are formulated in 
the last section. 



II. THE CARG METHOD 

To describe the method, we first introduce the nota- 
tion. We consider a model defined in a square lattice of 
lattice spacing a, and linear size L, with periodic bound- 
ary conditions (p.b.c). The spin variables ui are defined 
on each site i of the lattice and the Hamiltonian has the 
form 



H 



^ ' K a S a 



(1) 



where each S a is some kind of combination of the spin 
variables, with the only requirement of translation invari- 
ance subject to p.b.c. For our study we will consider the 
2D Ising model, including up to three even interactions 
(nearest-neighbor, second-neighbor, and four spin) and 
one odd interaction (a weak magnetic field) . 

The aim of RG theory is to study the critical properties 
of a model. The critical exponents for example, can be 
obtained from the linearized RG transformation matrix 
T a ii* , defined in eqn. (jHJ), by computing its eigenvalues. 
T a p* can be obtained numerically from the coupled equa- 
tions, 



Here, the super-indices (n— 1) and (n) denote the orig- 
inal and the renormalized quantity respectively, after one 
RG transformation, and the thermal average ( F(4>) ) is 
defined by the usual formula 



(3) 



The left and right quantities appearing in eqn. ([2]) can 
be obtained through the identities 



-( S 7 (n) Sp^-V ) + ( S* 7 (,i) ) ( Sp^-V ) 

(4) 



= -{ S* 7 (n) Sj n) ) + ( S 7 (n) ) ( Sj n) } (5) 



Due to the definition of the renormalized Halmiltonian, 
the thermal expectation of any function F of the variables 
a ( n ) yield the same value, whether one evaluates it using 
or fft"- 1 ) i.e.: 



Z(") 



(6) 
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Therefore, the quantities appearing in equations (j4|) 
and ([5]) can be computed directly by a numerical sim- 
ulation of the original Hamiltonian on the fundamental 
lattice. Because the simulations must be performed close 
to the critical point of the system, the configurations pro- 
duced by a local-update algorithm are not statistically in- 
dependent, which leads to inaccurate thermal averages. 
This phenomenon is called in the literature critical slow- 
ing down Q. In order to avoid this effect we use the 
collective-mode algorithm developed by Wolff @ for the 
Ising model. 

The expectation values appearing in equations and 
([5]) should be computed by using a simulation with a 
cluster algorithm on the fundamental lattice. These nu- 
merical values inserted into eqn. |2]) lead to coupled alge- 
braic equations, which allow to obtain a sequence of ap- 
proximations to the linearized RG transformation T a p*, 
as the renormalization transformation is iterated to the 
fixed point Hamiltonian H* , defined in the vicinity of the 
fixed point by 



K, 



(n+l) 



K P * ) (7) 



where the linearized RG transformation is defined by 



Ti * 
aP — 



dKp^- 1 ) 



(8) 



To evaluate T a p* one in principle needs a "linear re- 
gion" close to the critical temperature, where its deriva- 
tives are essentially constant. This region can be found 
by using FSS analysis of the lattice shifted critical tem- 
perature Q, and in the case of the 2D Ising model, by 
using the Binder Cumulant (see section IV). 

Finally, the critical exponents are obtained from the 
eigenvalues of T a p* from eqn. ((8]) in the standard way. 
For the 2D Ising model they are given by the fundamental 
relations v = lns/lnAi e and rj = d + 2 — 21nAi°/ln6, 
where Ai 6 ' - 1 is the largest even (odd) eigenvalue. 
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III. UNIVERSAL FLUCTUATIONS IN THE 2D 
ISING MODEL 

The 2 dimensional Ising model is well known and it is 
defined by the Hamiltonian 

<i,j> 

where Oi is the spin variable defined on the lattice site 
i, J > is the ferromagnetic constant, ks is the Boltz- 
mann's constant and 53<i •> stands for sum over nearest 
neighbors < i,j >. We use a system of units where 
Boltzmann's constant is set equal to unity throughout 
the paper and identify T with the reduced temperature 
T/J. It is well known that for an infinite square lattice 
(L — > oo), this model has a second order phase transition 
at the critical temperature T c defined by sinh(l/T c ) = 1 
or equivalently T c — 2/^n(l + y / (2)), which was first com- 
puted by Onsager [l7l |. 

Universality of Rare Fluctuations in Turbulence and 
Critical Phenomena, as it was first proposed in [l[ , linked 
two quite different physical systems: a confined turbulent 
flow and the finite volume 2D XY model in the low tem- 
perature regime. The (generalized) universality of both 
systems was based on the collapse of the probability dis- 
tribution functions of the corresponding global quanti- 
ties, the power consumption in the turbulent experiment 
and the magnetization in the critical system. The PDFs 
were conveniently normalized to their first two moments, 
and felt onto a common curve, which was called "univer- 
sal fluctuations" . 

A further generalization of the phenomenon described 
above was performed in [§], including non-equilibrium 
models like the auto igniting forest fire model and the 
Bak-Tang-Wiesenfeld sandpile model, and the equilib- 
rium 2D Ising model at temperature T*(L), which they 
called "bulk temperature" . 

Nevertheless, using the 2D and 3D Ising model Zheng 
et al. Q showed a dependence of the "universal fluctu- 
ations" on the equivalence class of the model. In a high 
precision MC simulation of the finite volume 2D XY- 
model, a slight but systematic dependence of the "uni- 
versal fluctuations" on the system temperature was first 
suggested [l(| • This claim was proved analytically in [12] • 
Moreover, by computing the skewness and the curtosis in 
the harmonic 2D XY-model (the second and third mo- 
ments of the PDF) the temperature dependence of the 
PDF was confirmed numerically and analytically in refs. 
fl3l | and [3] respectively. 

The precise physical definition of T*(L) is still in 
progress. In more recent papers, T*(L) has been linked 
to intermittency of the magnetization [151 ] . 

To address this issue we used the CARG method de- 
scribed in the above section and we have simulated the 
2D Ising model in a square lattice of lattice sizes L = 45, 
L — 64 and L = 108, with periodic boundary conditions, 



in a range of temperatures within the bulk T*(L) and 
the lattice shifted critical temperature T C (L) which can 
be defined, for example, as the temperature at which the 
magnetic susceptibility has a peak. 
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FIG. 1. The Binder Cumulant is shown for different lattice 
sizes. From their interception the critical temperature for 
infinite volume, T c (oo) = 2.26903 ± 0.00059, is obtained. It 
differs in 0.7 per one thousand from Onsager's exact result. 



IV. NUMERICAL RESULTS 

As it was mentioned in section II, and estimation for 
the critical temperature for the 2D Isin g m odel can be 
obtained by using the Binder Cumulant [181 ]: 



^) = !-3(W (10) 

where (M) is the averaged magnetization. It follows 
that u(L, T) depends on the system size and temperature 
through the ratio of the averaged fourth power to the av- 
eraged square of the magnetization. Nevertheless, from 
finite size scaling analysis, close to the critical tempera- 
ture it is independent of the lattice size, and therefore, 
the curves representing the Binder cumulant correspond- 
ing to different lattice sizes must collapse onto one point. 
This behavior is displayed in Fig. [TJ The value obtained 
agrees with Onsager's critical temperature within an er- 
ror of 0.7 per one thousand. 

Now we use CARG method to compute the critical 
exponents v and 77, associated to the Ising model. In or- 
der to compare with the results obtained in we use 
a lattice of lattice size L = 45 to compute the quanti- 
ties appearing in the RG equations (0])-([5]) close to the 
critical temperature, obtained as explained above, and 
perform two RG steps. Using the Wolff algorithm [||, 
10 5 sweeps were used to thermalize the system and 10 6 
configurations were used to compute thermal averages. 

For the RG analysis up to three even interactions 
(nearest-neighbor, second-neighbor and four spin or pla- 
quette) and one odd interaction (magnetic field) were 
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considered. In order to compare our results with the cor- 
responding ones obtained by [1], we display in table Q] the 
results of a lattice of lattice size L = 45, where the scale 
factor 6 = 3 was used. After two RG steps, the values 
compare quite well. Compared to their exact values of 
the infinite volume limit, they agree up to 0.25% for v 
and 0.52% for rj. 

TABLE I. The critical exponents obtained by the CARG 
method [PZ] and by Swendsen's method [S] for L — 45 are 
displayed. 



RG Step 


Ai e 


Ai° 


V 


V 


Ms) 


2.887 


7.712 


1.036 


0.2812 


1[PZ] 


2.8900 


7.7408 


1.0352 


0.2744 


2[s] 


3.006 


7.835 


0.998 


0.2524 


2[pz] 


3.0083 


7.8508 


0.9975 


0.2487 


Exact 


3 


7.8452 


1 


0.250 



In table [IT] we display the results for the critical expo- 
nents v and T), for a lattice of lattice size L = 64. The 
scale factor 6=2 was used. The results are even bet- 
ter than with the smaller lattice, with errors of order less 
than 0.13%. The reason for this behavior is physically in- 
tuitive: An important feature of the RG method is that 
the smallest system considered should be still large com- 
pared to the range of the fixed-point Hamiltonian, so that 
any significant truncation should be avoided. The valid- 
ity of this requirement is improved with the larger lattice 
sizes. 

TABLE II. Critical exponents for L = 64 and T = 2.259. 



RG Step 


Ai e 


Ai° 


V 


V 


1 


1.9586 


3.6856 


1.0311 


0.2362 


2 


1.9986 


3.6676 


1.0010 


0.2503 


Exact 


2 


3.6680 


1 


0.250 



In spite of the good accuracy obtained for the critical 
exponents performing only two RG steps, we iterated the 
RG transformations to three RG steps in order to com- 
pare our results with the ones obtained in reference [l9l |. 
The comparison is displayed in table III. Both methods 
achieve a remarkable accuracy, compared to the exact 
results. The main difference is related to the computa- 
tional time needed to obtained these values. Indeed, it 
turned out that the CARG method is faster than Swend- 
sen's method by a factor which increases monotonically 
with the lattice size, from six for L — 64 until 10 for 
L = 162. 

We want now to address the question whether the 
model displays or not a critical behavior at the bulk tem- 
perature T*(L), which would lead to scale invariance of 



TABLE III. The critical exponents obtained by CARG 
method [PZ] and those obtained by Swendsen [lj| for L — 
108. 



Kb Step 


Ai 


Ai 


V 


V 


ks] 


2.852 


7.705 


1.048 


0.2828 


l\PZ] 


2.8635 


7.7062 


1.0443 


0.2825 


2[S] 


3.021 


7.828 


0.994 


0.2540 


2[FZ] 


3.0027 


7.8269 


0.9992 


0.2542 


3[s] 


3.007 


7.831 


0.998 


0.2534 


3[FZ] 


3.0013 


7.8361 


0.9996 


0.2521 


Exact 


3 


7.8452 


1 


0.250 



the system. By definition, at T*(L) the PDF of the mag- 
netization has a similar form to the distribution originally 
obtained for the magnetization in the 2D XY-model in 
the zero temperature limit [T3 | . The scale invariance is 
usually expressed as a power law behavior of the physi- 
cal quantities involved and therefore we use the CARG 
method to look for critical exponents in the whole inertial 
range [T*(L), TJL)} @. 

In reference |15|, T*(L) was defined as the temperature 
for which the skewness ( the third normalized moment of 
the PDF of the magnetization ) is equal to that for the 2D 
XY model in the low temperature region. According to 
this definition, and for a lattice of lattice size L — 64 the 
numerical value T*(L) = 2.11 was found. In this article it 
was argued that at T*(L) there are (strong) correlations 
on all scales up to a length of the order of the system 
size L, which is an important feature of critical behavior. 
Furthermore, the fluctuations of the magnetization were 
linked to intermittency. To address this issue we have 
simulated a lattice with L = 64 as well. 

In Fig. [2j the relative errors between the exact values 
of the infinite volume critical exponents v and 77 and the 
corresponding numerical result are displayed as a func- 
tion of the system temperature. They are expressed in 
percents and from this figure we conclude that the rel- 
ative errors increase monotonically as the temperature 
moves away from T C (L), which includes the particular 
value T*(L), at which the relative errors are large. It is 
therefore rather unjustified to expect scale invariance of 
the physical quantities close to T*(L). 

Now we study a possible scaling behavior of the order 
parameter -the magnetization- in the vicinity of both the 
bulk and the shifted critical temperatures (T*(L) and 
T c (L))respectively. It is well known that the magnetiza- 
tion near the critical temperature behaves as a power law 
of the kind 

(M) ~ tP (11) 
where the reduced temperature r is given by 

T J T ~^)\ (12 ) 
T C (L) V ' 
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FIG. 2. Relative errors for the critical exponents r\ and v ob- 
tained in the second renormalization step with lattice size 
L — 64 and NAV = 10 6 . The minimum error occur at 
T c (64) = 2.259 and the critical exponents for this temper- 
ature are shown in table HT1 



malized trajectory, see second graphic in fig. [3l For the 
present computation we used Swendsen- Wang's cluster 
algorithm Q to update the fundamental lattice. We used 
the block spin parameter b — 2 and T C (L = 64) = 2.3008 
for the lattice shifted critical temperature. 




and (5 is the critical exponent. A direct numerical com- 
putation of this exponent is rather subtle because of the 
reflection symmetry of the Hamiltonian © , which leads 
in numerical simulations close to the phase transition to 
the use of the absolute value of the magnetization instead 
of the magnetization itself. One can instead compute j3 
by using the scaling relation /3 = v{d + rj — 2)/2. 

In ref. [HI the stochastic evolution of the magnetiza- 
tion was computed and it was argued that the system dis- 
plays "larger fluctuations" at T*(L) than at T C (L). This 
apparent phenomenon has been called "intermittency" . 
Nevertheless, and according to the discussion associated 
to fig. [2] no underlying scale invariance is responsible 
for this phenomenon. This apparent behavior is rather 
related to the fact that, in order to avoid metastable 
states in numerical simulations, the observable used in 
MC simulations is the absolute value of the magnetiza- 
tion instead of the magnetization itself. This quantity is 
bounded from below and therefore at T C (L) its fluctua- 
tions appear suppressed by a factor of order 2 compared 
to the corresponding at T*(L), where this lower bound 
plays no role as the magnetization does not vanishes in 
its neighborhood. 



FIG. 3. PDFs (raw and normalized data) for the fundamental 
lattice (L = 64) at temperature T C (L) = 2.3008 and for the 
two coarse grained lattices (L = 32 and L = 16). 

The invariance of the PDF under RG transformations 
can be understood by observing that it can be written 
as a Fourier transform of a partition function of an aux- 
iliary theory, which differs from the original theory by a 
dimension zero perturbation, with a very small imaginary 
coefficient. This was first pointed out in [jjj]. From the 
RG theory we know that partition functions are invariant 
under RG transformations. 

Finally we plot in figure [4] the PDF in the whole "iner- 
tial range" of temperatures. The curves depend slightly 
on the temperature close to T*(L), but this particular 
value, as we have argued, is not related to scale invariance 
or power law behavior of the system. Moreover, different 
curves normalized to the first two moments are difficult 
to distinguish from one another,close to the maximum. 
This has been pointed out in [la ], where the PDF at 
the particular value T*(L) was fitted by a Gumbel dis- 
tribution, which plays an important role in Statistics of 
Extremes. 



RG-flow of the Probability Density Function 

In this subsection, we proceed to study of the PDF by 
using the CARG method. In particular we have com- 
puted the RG flow of the PDF starting on a lattice of 
lattice size L = 64, which we called fundamental lattice, 
which we denote by A. We have performed two RG trans- 
formations. Two further lattices were defined with lattice 
sizes L — 32 and L = 16, which are denoted by A' and A" 
respectively. In fig. [3] the PDF is displayed for the three 
lattices. It follows that, properly normalized to the first 
two moments, the PDF remains invariant along a renor- 



V. CONCLUSIONS 

In this paper we have introduced a very accurate 
method, which we call CARG method, to study criti- 
cal systems close to a critical point. Compared to the 
original method proposed by Swendsen, our method is 
faster by a factor which grows linearly with the lattice 
size from 6 for L = 64 until 10 for L = 162, due to the use 
of cluster algorithms to simulate the fundamental Hamil- 
tonian. This advantage allows to simulate larger lattices 
in reasonable computational times, and to improve the 
accuracy of the results, as shown in table III. 
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FIG. 4. The PDF is displayed for temperatures in the "inertial 
range" and for the lattice size L = 64. 

We illustrated the method by using the 2D Ising model 
defined in a square lattice of lattice size L, and we have 
shown how to obtain very accurate values for the critical 
exponents without the previous knowledge of the criti- 
cal temperature. We used further the 2D Ising model to 
study scale invariance and universality as the underlying 
mechanism which could give rise to the approximate gen- 



eralized universal behavior of fluctuations. We computed 
the probability density function (PDF) of the magneti- 
zation and its RG trajectory close to the lattice shifted 
critical temperature. 

Critical behavior is associated to scale invariance, 
which is commonly represented as a power law behav- 
ior: the physical quantities are described close enough to 
a critical point, by critical exponents, which characterize 
the equivalence class the system belongs to. In the spe- 
cial case of the 2D Ising model, from Fig. [5] we conclude 
that, at the bulk temperature T*(L), neither scale invari- 
ance nor universal behavior is really present. Therefore 
the approximate collapse of the PDF onto the "universal 
distribution" seems not to be related to critical behav- 
ior but rather to a numerical phenomenon associated to 
approximated scaling relations and extreme statistics, as 
proposed in [l6| , and to the constraint character of the 
global quantity used to compute the PDF. 
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